clear all;clc;close all;
cs=CubeSolver('nMax',11,'lambda',1.4,'nGrid',100,'lengthD',1.2,...
              'diffuseD',1.1,'fS',0.1,'fF',0.9,...
              'OD',5.2,'power',10.2,'P0',2.3);
w=0.2;
beamshape=@(x,y,z) exp(-((x-0.5).*(x-0.5) + (y-0.5).*(y-0.5))/2/w^2);%.*boolean(x>0.5).*boolean(y>0.5);
cs.SolveMz(beamshape,'cutOff',4);

cs_new=CubeSolver('nMax',11,'lambda',0.4,'nGrid',100,'lengthD',1.2,...
              'diffuseD',1.1,'fS',0.1,'fF',0.9,...
              'OD',5.2,'power',10.2,'P0',2.3);
vecVal=cs.Results.coeffMz;
cs.SolveT2(beamshape, vecVal,'cutOff',4,'tList',(0:0.01:0.5)');
%%
close all;
plot(cs.Results.tList,reshape(cs.Results.Mp(50,50,50,:),[],1));
xlabel('t/\tau_D');ylabel('Mp(0.5,0.5,0.5,t)');
cs.dispPolarize(beamshape);
cs.plotCubic(cs.Results.Mp(:,:,:,1));
cs.plotCubic(cs.Results.Mp(:,:,:,end));
